A novel defined pyroptosis-related gene signature predicts prognosis and correlates with the tumour immune microenvironment in lung adenocarcinoma

Lung adenocarcinoma (LUAD) is one of the most common causes of cancer-related death. The role of pyroptosis in LUAD remains unclear. Our study aimed to identify a prognostic signature of pyroptosis-related genes (PRGs) and explore the connection of PRGs with the tumour microenvironment in LUAD. Gene expression and clinical information were obtained from The Cancer Genome Atlas database. Consensus clustering was applied to classify LUAD patients. The least absolute shrinkage and selection operator Cox and multivariate Cox regression models were used to generate a PRG-related prognostic signature. The correlations between PRGs and tumour-infiltrating immune cells or the tumour mutational burden were analysed by Spearman’s correlation analysis. In this study, 44 PRGs significantly differed in expression between LUAD and normal tissues. Based on these genes, patients were clustered into three clusters with significantly different distributions of tumour-infiltrating immune cells and immune checkpoint regulators. A total of four PRGs (NLRP1, HMGB1, CYCS, and BAK1) were used to construct a prognostic model. Significant correlations were observed between these prognostic PRGs and immune cell infiltration or the tumour mutational burden. Predictive nomogram results showed that BAK1 could be an independent prognostic biomarker in LUAD. Additionally, the expression level of BAK1 was validated in two independent Gene Expression Omnibus cohorts. Our identified prognostic PRG signature may provide insight for future studies targeting pyroptosis and the tumour microenvironment in LUAD. Future studies are needed to verify our current findings.

The relationship between pyroptosis and cancer is complex since the release of inflammatory mediators can form a microenvironment suitable for tumour cell growth, while on the other hand, as a type of cell death, the induction of tumour pyroptosis can inhibit the occurrence and development of tumours 9,[12][13][14] . Therefore, pyroptosis is thought to play a dual role in tumours. In NSCLC, a higher expression level of GSDMD is related to invasive features, including an advanced tumour-node-metastasis stage and enlarged tumour size 12,15 . Downregulation of gasdermin A (GSDMA) results in caspase-3 activation and cancer cell death through the mitochondrial apoptotic pathway 16,17 . Researchers also found that GFNA5/GSDMD could determine the caspase-3-activated cell death mode and drug reactivity 12,18 . Some studies have explored the roles of pyroptosis-related genes (PRGs) in LUAD. Song et al. 19 constructed a pyroptosis-related lncRNA signature to predict prognosis in LUAD patients. Lin et al. 20 and Liu et al. 21 each applied the same pyroptotic gene set consisting of 33 PRGs to investigate the roles of these genes in LUAD. These studies provide convincing evidence that pyroptosis is closely connected with LUAD pathogenesis.
In this study, based on an improved pyroptosis-related gene set consisting of 52 genes, we clustered LUAD patients from The Cancer Genome Atlas (TCGA) database and explored the tumour immune status within each subclass. Then, we developed a prognostic signature using the least absolute shrinkage and selection operator (LASSO) Cox method and studied the correlations between prognostic PRGs and the tumour mutational burden (TMB) or tumour-infiltrating immune cells. The expression level of BAK1, a potential independent predictor in LUAD, was validated in two independent Gene Expression Omnibus (GEO) datasets. Our findings indicate the potential connections of pyroptosis with disease prognosis and the immune microenvironment in LUAD patients.

Materials and methods
Acquisition of gene expression and clinical data. We obtained the RNA-sequencing data of 513 LUAD patients and corresponding clinical information from TCGA database (https:// portal. gdc. cancer. gov/ repos itory). Data and figures were analyzed and generated by R software (version v4.0.3; https:// www.r-proje ct. org/).
Estimation of immune cell type enrichment. Immune cell enrichment analysis of the RNA-seq data was performed with xCell. Relative cell type abundance was quantified and visualized across all samples. The abundance of each cell type across the clusters was compared using the Kruskal-Wallis test. Cell types with p < 0.05 were considered significantly differentially enriched.
Development of a PRG-based prognostic model. Cox regression analysis was performed to evaluate the prognostic roles of PRGs. Kaplan-Meier survival analysis was applied to compare survival difference between two groups, with the p-values and hazard ratios (HRs) of 95% confidence intervals (CIs) generated by log-rank tests and univariate Cox proportional hazards regression. PRGs with significant prognostic values were included in further analyses. LASSO Cox regression analysis was used to construct a prognostic model based on the identified PRGs. The LUAD patients were divided into low-and high-risk clusters according to the median risk score, and the overall survival (OS) time was compared by Kaplan-Meier analysis. Predictive accuracy was evaluated by performing time receiver operating characteristic (ROC) curve analysis.

Construction of a gene-based prognostic nomogram. A composite nomogram was constructed
based on the results of multivariate Cox proportional hazards analysis to predict 1-, 3-, and 5-year overall recurrence. A forest plot was used to present the p-value, HR, and 95% CI of each variable via the 'forestplot' R package.

Immune infiltration and tumour mutational burden (TMB) analysis. The correlation between each
prognostic PRG and tumour-infiltrating immune cells was analysed by Tumour Immune Estimation Resource (TIMER, https:// cistr ome. shiny apps. io/ timer/), a web portal for comprehensive analysis of tumour-infiltrating immune cells. Spearman's correlation analysis was applied to describe the correlation between each prognostic PRG and the TMB. A p-value less than 0.05 was considered statistically significant.

Results
Identification of differentially expressed PRGs between normal and LUAD tissues. The expression levels of 53 PRGs were compared between 49 normal and 513 LUAD tissues from the TCGA. Fortyfour PRGs were significantly differentially expressed (Fig. 1A)

Identification of LUAD clusters by consensus clustering. To explore the connections between the 44
PRGs and LUAD subtypes, we performed consensus clustering analysis of all 513 LUAD patients in the TCGA cohort. The number of clusters was represented by the parameter 'k' . The empirical cumulative distribution function (CDF) was plotted to determine the optimum k value for the sample distribution to reach maximal stability ( Fig. 2A). By increasing k from 2 to 6, we found that when k = 2, LUAD patients could be divided into three distinct and nonoverlapping clusters (Fig. 2B). The differential gene expression profile is presented in a heatmap (Fig. 2C). The OS time was compared among the three clusters, and no obvious differences were found (p = 0.39, Fig. 2D).
Distinct TMEs of PRG-based clusters. Next, we analysed tumour-infiltrating immune cell data and found that most of the infiltrating immune cells were abundant in Cluster 3, which included T cells (naive CD8 + T cells, effector memory CD8 + T cells, central memory CD8 + T cells, central memory CD4 + T cells, and effector memory CD4 + T cells), B cells (naive B cells, plasma B cells, class-switched memory B cells, and memory B cells), macrophages (M1 and M2), dendritic cells (activated myeloid dendritic cells and plasmacytoid dendritic cells), monocytes, mast cells, neutrophils and eosinophils (Fig. 3A). We also found that the proportions of immune cells within each of the three clusters were different, while the composition types were the same (Fig. 3B). We also checked the expression of immune checkpoint-related genes (SIGLEC15, CD274, HAVCR2, PDCD1LG2, CTLA4, TIGIT, LAG3, and PDCD1) within the 3 clusters. The results revealed an expression profile similar to that of tumour-infiltrating immune cells in the 3 clusters. As shown in the immune checkpoint-related gene expression heatmap, where different colours represent the expression trends in each sample, the expression of immune checkpoint-related genes was most enriched in Cluster 3 (Fig. 3C).

Construction of a PRG-based prognostic model.
To construct a PRG-based prognostic model, univariate Cox expression analysis was performed to screen the differentially expressed PRGs for prognostic value. We identified 8 genes with significant prognostic value in Kaplan-Meier survival curves (Fig. 4). The results suggested that a poor survival rate in LUAD patients was related to high expression of BAK1 ( with the most significant p-values (Fig. 5A,B). The risk score = (0.1918) *BAK1 + (0.2309) *CYCS + (0.0635) * HMGB1 + (− 0.1105) * NLRP1. Based on the risk score, LUAD patients were separated into two clusters. The risk score, overall survival time and expression of these four PRGs are presented in Fig. 5C. As the risk score increased, the survival time of patients decreased. HMGB1, CYCS, and BAK1 were associated with an increased risk score, while NLRP1 was a protective gene that was correlated with a decreased risk score (Fig. 5C). Kaplan-Meier curves indicated that patients with a high risk score had a worse overall survival probability than those with a low risk score (median time: 3.2 years in the high-risk subgroup versus 7.1 years in the low-risk subgroup, p < 0.001, Fig. 5D). Time-dependent receiver operating characteristic (ROC) curve analysis was applied to evaluate the sensitivity and specificity of the prognostic model, and the results showed that the areas under the ROC curves (AUCs) were 0.646 for 1-year survival, 0.592 for 2-year survival and 0.662 for 5-year survival (Fig. 5E). we built a nomogram including the clinical features (age, sex, race, pTNM stage, and smoking status) that were generally believed to have a certain impact on the prognosis of LUAD and the 4 prognostic PRGs to predict the survival rate of LUAD patients (Fig. 6). Univariate and multivariate analyses indicated that BAK1 expression might be a candidate independent factor, similar to pTNM stage, that affects the prognosis of LUAD patients (Fig. 6A,B). A C-index of 0.693 indicated that the nomogram had good predictive value (Fig. 6C). The predictive nomogram suggested that 2-year, 3-year and 5-year overall survival rates could be predicted relatively well according to an ideal model in the entire cohort (Fig. 6D).

Associations of the tumour mutational burden (TMB) and tumour-infiltrating immune cells with prognostic PRGs.
The TMB has been identified as a promising biomarker for predicting immunotherapy responses in patients with NSCLC. To explore the possibility of using prognostic PRGs as biomarkers in immunotherapy for LUAD, we analysed the correlations of prognostic PRGs with the TMB (Fig. 7). The results revealed positive correlations between the TMB and CYCS (Fig. 7B, p = 1.54e −5 ) or HMGB1 (Fig. 7C, p = 0.04) and a negative correlation between TMB and NLRP1 (Fig. 7D, p = 2.1e −6 ). However, there was no significant correlation between TMB and BAK1 (Fig. 7A, p = 0.96). Furthermore, to determine the roles of these prognostic PRGs in the tumour immune microenvironment, we performed correlation analysis of tumour-infiltrating immune cells with each prognostic PRG in LUAD by the TIMER database (Fig. 7E). Our results showed positive correlations between NLRP1 expression and the abundance of B cells (p < 0.01), macrophages (p < 0.01), myeloid DCs (p < 0.01), neutrophils (p < 0.01), CD4 + T cells (p < 0.01) and CD8 + T cells (p < 0.01). Significant negative correlations were observed between the expression of HMGB1 and B cells (p < 0.01), myeloid DCs (p < 0.01), and CD4 + T cells (p < 0.01). The same negative correlation www.nature.com/scientificreports/ pattern was also observed for CYCS with immune cells, including B cells (p < 0.01), myeloid DCs (p < 0.01), CD4 + T cells (p < 0.01), macrophages (p < 0.01) and neutrophils (p < 0.05). We also found positive correlations between BAK1 expression and myeloid DCs (p < 0.05) or neutrophils (p < 0.01).
Construction of a prognostic model with BAK1. As indicated by our analysis of the predictive nomogram results, we concluded that BAK1 might be an independent prognostic biomarker in LUAD. We further performed prognostic analysis of BAK1 in patients based on the risk score. The risk score distribution, survival time and BAK1 expression are shown in Fig. 8A. As the expression of BAK1 increased, the risk score  www.nature.com/scientificreports/ increased; thus, patients' risk of death increased accordingly (Fig. 8A). Kaplan-Meier curves revealed that patients with a high risk score presented a worse overall survival probability than those with a low risk score ( Fig. 8B; median time: high-risk group, 3.5 versus low-risk group, 4.4; p = 0.0137). The AUCs of the risk score were 0.552 for 1 year, 0.543 for 3 years, 0.619 for 5 years and 0.763 for 10 years (Fig. 8C).

Validation of BAK1 expression in two independent GEO datasets. Two independent GEO cohorts
were utilized to validate the BAK1 expression results in LUAD patients. GSE 10,799 included 3 normal controls and 16 LUAD patients, while GSE 66,759 included 5 controls and 76 LUAD patients. The Wilcox test was used to compare BAK1 expression between LUAD and normal samples. The results consistently showed that BAK1 expression was significantly higher in LUAD samples than in control samples (Fig. 9).

Discussion
Pyroptosis is a newly identified type of programmed cell death that plays a dual role in tumour development and therapeutic mechanisms 7,27 . The process is characterized by rapid rupture of the cell membrane and release of proinflammatory intracellular contents 8,9,13,14 . This unique type of death has led to considerable studies on pyroptosis in various tumours [28][29][30][31] . However, the specific role of pyroptosis in LUAD remains unclear.
In this study, we explored 52 currently known PRGs in LUAD and normal tissues and identified that most of them (44/53) were differentially expressed. The three clusters defined based on the differential PRGs did not show any significant differences in overall survival time. However, the abundance of tumour-infiltrating immune cells and immune checkpoint regulators showed significant differences among the three clusters. To further assess the www.nature.com/scientificreports/ prognostic value of these PRGs, we constructed a 4-gene risk signature via Cox univariate analysis and LASSO Cox regression analysis and found that it had good accuracy for predicting the survival of LUAD patients. Survival analysis and prognostic models were developed based on four genes (BAK1, CYCS, HMGB1, and NLRP1) in this study. Among these, BAK1 was an independent risk factor for OS in patients with LUAD. High levels of endogenous BAK have been observed in both small cell lung cancer and NSCLC cell lines. In addition, increased BAK expression was correlated with a poor prognosis in NSCLC patients 32 . Recent research has suggested that BAK could be a promising prognostic indicator and potential therapeutic target in lung cancer patients 33 . These studies are in line with our findings, indicating that BAK1 has a positive correlation with LUAD patient prognosis. We confirmed the high BAK1 expression level in two independent GEO cohorts. Pyroptosis in colon cancer cell lines can be mediated by BAK1 or BAX alone, and caspase 3 activity is required in BAK/ BAX-mediated pyroptosis 34 . It is reasonable to hypothesize that abnormal BAK1 expression could promote the development of LUAD by regulating the pyroptotic pathway 34 . However, the mechanisms by which BAK1 mediates pyroptosis and leads to LUAD remain elusive. CYCS, as a mitochondrial protein that participates in the regulation of cell death, was reported to be an oncogene in LUAD 35 . Previous research has demonstrated that serum CYCS levels are correlated with disease progression and a poor prognosis in NSCLC 36 , which is consistent with our current results. Feng et al. proved that HMGB1 was overexpressed in NSCLC tissues 37 . We also found that higher expression of HMGB1 was correlated with a poor clinical prognosis in LUAD patients. Similar results from Chang et al. showed that HMGB/RAGE signalling was significantly associated with patient prognosis, which agrees with our survival analysis results for the TCGA datasets 38 . NLRP1 is recognized as an important component of complexes that can activate caspase-1 directly 39 . Shen et al. reported that the NLRP1 expression in LUAD tissue was considerably lower than that in normal tissues 40 . This decreased NLRP1 expression was associated with high T and N stages 40 . Consistently, we found that LUAD patients with low NLRP1 expression had a worse prognosis than those with high expression. Caspase-6 has been shown to be an important regulator in inflammasome activation that could promote the activation of programmed cell death pathways, including pyroptosis, apoptosis and necroptosis (PANoptosis) 41 . Studies have indicated that NLRP3 inflammasome activation enhances the proliferation and metastasis of the lung adenocarcinoma cell line A549, which are mediated by AKT, ERK1/2, and CREB, and upregulation of SNAIL 42 . Consistent with these findings, in our current study, the expression of CASP6 was negatively correlated with patient OS.
Recent research has demonstrated that the TMB is a determinant of immune-related survival in a variety of tumours, such as breast cancer and lung cancer [43][44][45][46] . In the current study, the prognostic pyroptosis-related genes www.nature.com/scientificreports/ CYCS and HMGB1 were positively correlated with the TMB, while the expression of NLRP1 was negatively correlated with the TMB. Pyroptosis was initially found in macrophages but has recently been identified in a variety of immune cells 8,11,47 . Another interesting finding in our results is that the above four prognostic PRGs (BAK1, CYCS, HMGB1, and NLRP1) were significantly correlated with immune cell infiltration, which further confirmed that pyroptosis in immune cells might participate in regulating the tumour microenvironment. Our study has great clinical significance, especially the prognostic model that can be used to predict the prognosis of LUAD patients. Although various advances have been shown to be beneficial to some patients, such as immune checkpoint therapies including programmed cell death protein-1 (PD-1) and T-lymphocyte-associated antigen 4 (CTLA4) blockade, a proportion of patients are resistant to current therapeutic strategies, partly due to the heterogeneity in PD-L1 expression, the TMB and T-cell infiltration in LUAD 48 . Clinically, PRG-based classifiers have the potential to provide a novel approach for identifying novel subtypes of LUAD and personalized treatment for these patients. In terms of immune checkpoint therapy, some of these agents could affect the immune microenvironment through pyroptosis. Therefore, TMB-related PRGs have the potential to be used to guide the curative efficacy of immune checkpoint inhibitors. This study provides new insight into the molecular mechanism underlying LUAD pathogenesis. www.nature.com/scientificreports/ Our current study is an exploratory analysis conducted using the TCGA-LUAD cohort, and the results presented here will need to be confirmed in larger datasets. Additionally, the results presented here will need to be confirmed by in vivo and in vitro experiments.
In conclusion, our study identified 44 PRGs that were differentially expressed between LUAD and normal tissues. Based on these PRGs, LUAD patients were classified into 3 subclusters with differential expression levels of tumour-infiltrating immune cells and immune checkpoint regulators. A novel prognostic model based on four PRGs was constructed and used to predict the prognosis of LUAD patients. The correlations of PRGs with the TMB and immune cells may provide evidence that pyroptosis might play an important role in the tumour microenvironment.